Computation of 2-D Spectra Assisted by Compressed Sampling 
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The computation of scientific data can be very time consuming even if they are ultimately deter- 
mined by a small number of parameters. The principle of compressed sampling suggests that for 
typical data we can achieve a considerable decrease in the computation time by avoiding the need 
to sample the full data set. We demonstrate the usefulness of this approach at the hand of 2-D 
spectra in the context of ultra-fast non-linear spectroscopy of biological systems where numerical 
calculations are highly challenging due to the considerable computational effort involved in obtaining 
individual data points. 



Introduction - Signals of interest that are obtained in 
experiments or by numerical computation of natural phe- 
nomena are not white noise but contain hidden structure 
and redundancy. The fact that it is possible to make 
use of these structures is the basis of compressed sam- 
pling [IrH] which demonstrates, surprisingly, that the 
accurate reconstruction of such signals can be achieved 
with high probability from a number of data points that 
would be deemed insufficient by the Nyquist - Shannon 
criterion. This promise of considerable efficiency gains 
together with the development of numerical algorithms 
that allow for its efficient implementation have led to an 
explosion of applications of compressed sampling in sig- 
nal processing [5]. 

Here we apply the principle of compressed sampling to 
the problem of non-linear two-dimensional (2-D) spec- 
troscopy and demonstrate that its application can offer 
a considerable reduction in the computational effort in- 
volved in the numerical determination of such spectra. 
This makes accessible 2-D spectra for larger systems or 
for systems that are governed by more complex dynam- 
ical equations, for example due to the inclusion of non- 
Markovian features [5H5]- This computational saving is 
not restricted to 2-D spectra. It serves however to high- 
light the observation that numerical computation of sig- 
nals arising in nature can be made considerably more ef- 
ficient employing the paradigm of compressed sampling. 
This is true whenever the signal is determined by a small 
number of parameters but the computation of individual 
data points is very costly [101 E] ■ It is noteworthy that 
the same ideas may also be applied to reduce the exper- 
imental effort in measuring 2-D spectra. Here however 
the impact of experimental noise on the quality of this 
reconstruction requires careful consideration |12j and will 
be studied elsewhere. 

Principles - In non-linear 2-D spectroscopy a series of 
three (or more) laser pulses is used to probe the dynamics 
of a system and the resulting information is arranged as 
a matrix representing the two dimensional Fourier trans- 
form of the signal which in turn describes the system re- 
sponse in dependence of the spacing between laser pulses 



(more details will be given in the next section). The 
point by point computation of the entries of this matrix 
can be very time consuming, especially when the rele- 
vant dynamical system is high-dimensional or subject to 
a complex non-Markovian system-environment interac- 
tion. As explained above, we expect however that this 
matrix contains structure. In line then with the paradigm 
of compressed sampling it will be possible to achieve a 
considerable reduction in the number of elements of the 
signal that need to be obtained through (random) sam- 
pling (by experimental observation or numerical compu- 
tation) . After an outline of the numerical algorithm that 
is being used we provide a demonstration of the achiev- 
able computational gain at the hand of a specific bio- 
molecular complex of current interest. 

Let M be an n x n matrix that represents the 2-D spec- 
trum that we wish to reconstruct and let experimental 
observation or numerical computation determine a small 
subset of these elements described by the set of indices O. 
Denoting by tr the trace norm of a matrix X it can 
be proven that the solution of the minimization problem 

min[tr|A| : X t] = My for £ fi] (1) 

is unique and yields the matrix M with a probability 
larger than 1 — exp(—f3) if the number of sampled entries 
of M is of order nr(\ + j3) Inn, where f3 > is an ar- 
bitrary constant and r is the rank of the matrix M (see 
[l3l [14] for proofs and a rigorous mathematical state- 
ment). While reconstruction of the matrix M cannot 
be guaranteed with certainty, the probability for obtain- 
ing the correct reconstruction of the matrix M can be 
increased arbitrarily by repeating the sampling and re- 
construction procedure k times. Accepting the majority 
result increases the probability that the so obtained re- 
sult is correct to above 1 — e~ k/3 . For typical values of n 
in the range of 100 — 1000 only /3 = Inn this reduces the 
probability of a wrong reconstruction to 10 -8 — 10 -12 . 
Given that the computation time for the full data set 
scales as n 2 this suggest that a considerable potential 
reduction in computational effort is possible. 

The solution of the minimization problem eq. (fTl) can 
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be obtained by standard solvers for semidefinite opti- 
mization problems which do however tend to be limited 
to relatively small matrix sizes with several hundred el- 
ements only. Here we solve the minimization problem of 
cq. ([lj by means of the so-called singular value thresh- 
olding (SVT) algorithm [TSUIS] which permits very large 
matrices to be treated and to provably approximate the 
solution of eq. ([I]) to arbitrary precision. The SVT- 
algorithm solves iteratively the set of equations 

X {k) = v {k-i) max (£,(fc-i) _ T%) )i/( fe - 1 ) (3) 
y(fc) = y(fc-i) + SV n (M - X {k) ) (4) 

where Vq(M) is a matrix whose entries satisfy 
(Vn(M))ij — Mij for <E fl and are zero otherwise 

while r and S are constants whose choice will be described 
below. The algorithm is initialized by a random matrix 
for the first iteration k = 1. Then the singular value 
decomposition of Y^ ' is computed in eq. (|]). This is 
followed by the so-called soft thresholding step eq. |3| 
where we subtract from all the singular values the param- 
eter t and set negative results to 0. This step yields a ma- 
trix that has lower rank than y( fc_1 ). The final step 
of iteration k then uses X( k > to construct a matrix Y^ 
that satisfies the constraints Xy = My for E ft 

more closely than the matrix yC 0-1 ) cq. Q. These steps 
are iterated until convergence has been achieved. 

The choice S < 2 ensure provable convergence of the 
iteration. For r the choice r = 5n for n x n-matrices leads 
to fast convergence and very close approximation to the 
solution of eq. (see [15] for details) . Making use of the 
sparseness of the matrices the SVT algorithm is capable 
of treating very large matrices (reaching 30000 x 30000 
and above). Variants of this algorithm can achieve a con- 
siderable additional increase in computational efficiency 
[17[ 118] by employing linear time algorithms for the im- 
plementation of the singular value decomposition in eq. 

Whenever the calculation of individual data points is 
costly, it will be advantageous to apply the compressed 
sampling paradigm and compute only a small subset of 
the matrix elements and complete the remaining entries 
of the matrix with the SVT algorithm. In the following 
we would like to demonstrate the power of this approach 
in the calculation of 2-D spectra. 

Applications - Non-linear spectroscopy has proven to 
be useful in order to unveil the dynamics involved in ex- 
citonic transfer of light harvesting complexes, due to the 
fact that it is sensitive to excitonic quantum superpo- 
sitions, i.e., excitonic coherences [2D ~23 as well as vi- 
brational features of the protein environment [24] [25]. 
Non-linear 2-D spectroscopy can resolve the third-order 
polarization of the electronic system, from the signal 
S^ 3 \ti, £2, £3), arising from photo-excitation of three con- 
secutive pulses with wave vectors ki, k2 and k3, sepa- 



rated by time intervals t\ and ti and measured at a time 
£3 after the last pulse. It is customary to Fourier trans- 
form the time t\ and t 3 dimensions to yield S (uji, t^u^) 
in order to generate a 2-D spectra parametrized by the 
waiting time ti- The numerical computation as well as 
the experimental determination of these 2-D spectra rep- 
resent a challenging task because of the large number of 
measurements that need to be taken and the considerable 
computational resources required for the determination 
of S^(t u t 2 ,t 3 ). 

Here we will demonstrate by means of numerical 
examples, that the paradigm of compressed sampling 
delivers a considerable gain in computational efficiency 
as it is not necessary to determine the full signal 
S^ 3 )^!,^)^) for a given choice of t^. The computation 
for a randomly chosen subset of times t\,t$ suffice to 
reconstruct the full signal. 

We consider the 2-D spectra of a pigment-protein com- 
plex, the Fenna-Matthews-Olson (FMO) complex, which 
has received considerable attention recently 23 . In a 
spectroscopy experiment the pigment-protein complex 
dynamics is probed by means of applying different laser 
pulses and measuring the response to them. Using Liou- 
ville space notation we can describe the evolution of our 
system under the effect of the laser pulses by the master 
equation 

Here the Liouvillian Jzf describes the internal excitation 
dynamics including dephasing of the complex while the 
term Jz?i n t(t) stands for the interaction with the laser. 
Typically in these complexes the dissipation time scales 
are much longer than the dephasing induced by the en- 
vironment and hence we will consider only pure dephas- 
ing contributions to the internal dynamics. That is, we 
will consider an unperturbed Liouvillian with the form 
5£ = Jz?Ham + ^cicph • Within the Born-Markov approx- 
imation we can write each of the terms in the equations 
above as: 

J2kam= (6) 
7 

J% cph ^^2 7i KVf-p) (7) 

^i n t(t) = [H iDt (t),p], (8) 

with o~ x , o~ y , <J Z the standard Pauli matrices and 
Hint — —E(t) ■ V. For the homogeneous dephasing 
rates we have chosen random values within a realis- 
tic interval for each chromophore, i.e., {7i,}i=i..7 = 
{0.94, 2.26, 2.83, 1.70, 1.70, 1.88, 2.07} x 10- 3 rad/fs (with 
the conversion lrad • fs^/rcnT 1 = 2tt ■ 2.99792458 x 
10~ 5 ), which correspond to dephasing times in the scale 
of few picoseconds. As regards the actual values of 
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the site energies and intersite coupling rates in the one- 



exciton sector we have used the following Hamiltonian[! 
(in units of rad/fs): 
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Finally, the dipolar operator is defined from the indi- 
vidual dipole momenta of each chromophore {fii\i=\...-! 
as V — 2l=i fccf- The directions of the seven induced- 
transition dipoles fli, extracted from the structure of the 
FMO complex [28], are 
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(10) 



and their magnitude is taken to be the same and given 
in units of the dipole strength |/iBCh| of an isolated bac- 
teriochlorophyll. 

The signal measured on a general spectroscopic experi- 
ment is intimately related to the system electric response 
functions of different orders S'"', In particular in non- 
linear 2-D spectroscopy experiments the signal can be re- 
lated to the third-order response function (t\,t2,t3). 
This function can be written in a very compact form us- 
ing tetradic notation in Liouville space |20j 



5 (3) (ii,i 2 ,i 3 



((v\&(t 3 )y<#(t 2 yy&(t 1 )r\p(-oo)}) 

(ii) 

where ^{t) is the Liouville space Green Function in the 
absence of the radiation field 



Sf(t) = 6{t) exp (^*) 



(12) 



and the action of the superoperator Y upon an ordinary 
operator A is defined as 



VA = [V,A]. 



(13) 



The dynamics of a photo-reactive molecule interact- 
ing with a series of laser beams does not preserve the 
number of excitations within the molecule. If we con- 
sider the natural assumption that the system is in its 
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ground state when the first pulse arrives at the sample 
(i.e., p{— oo) = («? | ) , then the only relevant sectors to 
compute the third-order response function are the ground 
state itself, the one-exciton and the two-excitons mani- 
folds. In FMO, with seven chromophores, this accounts 
for 1 + 7 + 21 = 29 states, which means that any op- 
erator in Hilbert space can be represented as a matrix 
of dimensions 29 x 29, while a superoperator acting on 
the Liouville space can be represented by matrices of di- 
mension 29 2 x 29 2 . On the other hand, the only op- 
erators needed to compute the response function from 
cq. ( 11 ) are the total Hamiltonian of the FMO molecule 
(and the corresponding Liouvillian Jz?kam), 



the dephas- 

ing superoperator Jz?dcph and the total dipolar operator V 
(with the corresponding superoperator Y). It is indeed 
a straightforward procedure to construct these operators 
in the manifolds mentioned above [2 6) and it should be 
noted that the only information required to construct 
them are the physical matrix elements given by H s in eq. 
( 9j)) together with the dipole information contained in eq. 
(10) and the dephasing rates {7i}i=i...7 already written 
above. 

Results of the application of compressive sampling on 
these 2-D spectra are presented in Fig. [T] for population 
time ti — 0. Already computing only about 0.16% ran- 
domly sampled points (topleft) the position of peaks start 
to emerge. For higher sampling ratios of 1% (bottomlcft) 
many features of the spectrum are reproduced qualita- 
tively even though some minor features are missing. For 
8% (topright) the spectra are reproduced very well quan- 
titatively and the differences with the exact 2-D spectra 
(bottomright) are almost negligible as quantified by the 
Frobenius norm difference (see caption of Fig. [IJ . Note 
that all the spectra have been normalized to the same 
maximum value for ease of comparison. The reconstruc- 
tion of the 600 x 600-matrices from the randomly sam- 
pled points was possible within less than 1 minute on 
a standard laptop for S ( - 3) (t ll t 2 = 0,t 3 ) [29 . The full 
computation of such a spectrum on the same machine 
takes around 60 minutes and scales as n 2 for annxn 
matrix. Assuming a sampling of 1% of the data set the 
same calculation could have been achieved in around 2 
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FIG. 1. The 2-D spectrum of the FMO complex with — (with the parameters described in the text) for a random sampling 
of (topleft) 0.16%, (bottomleft) 1%, (topright) 8% of the total data set S 1 ' 3 ' (ti, ti = 0, and the exact spectrum (bottomright). 
The black lines indicate the eigenfrequencies of the Hamiltonian in the first exciton manifold. Even for very small sampling rates 
the principal features of the spectrum are well reproduced and for a sampling of 8% only very minor differences remain. Each 
time variable has 600 bins resulting in a total data set in the form of a 600 x 600 matrix. The convergence of the reconstructed 
matrices Sk to the exact spectrum S EX with increasing sampling k can also be quantified using the standard Frobenius norm 
with the magnitude T> k = \\S C x — SfcH/HS'exll which yields {V k } k =o.i6, 1.00,8.00 = {1.81,0.38,0.08}. For the sake of comparison, 
the same measure averaged over a large sample of entirely random matrices S ia .nd yields £^ ran d = 8.96. On the other hand, the 
most distant matrix from the exact spectrum g ives i^worst — 15.72. 



minutes including the SVT reconstruction implying a 30 
fold improvement in computation time. The relative ef- 
ficiency improvement grows almost linearly with n. It 
should also be noted that 2-D spectra simulations aimed 
at direct comparison with experiments typically include 
rotational and inhomogeneous averages that would mul- 
tiply the total computation time by the size of the sam- 
pling distribution used to make each independent aver- 
aging. This does not affect our general conclusion con- 
cerning the computational speedup that can be obtained 
by means of compressive sampling. 

This example shows that even for rich 2-D spectra such 
as those obtained for the FMO complex, compressive 
sampling yields a considerable reduction of the number of 
arguments ti,t 3 for the which S^(ti,t 2 = 0,^) needs to 
be evaluated to obtain the complete 2-D spectrum with a 
precision that is sufficient for comparison to experiment. 
This exemplifies the power of signal processing concepts 
to save computational resources. The extraction of in- 



formation from complex spectra can also be assisted by 
methods from signal processing and in a future publica- 
tion we will explore the use of wavelets to this effect. 

Conclusions - We have demonstrated that the princi- 
ples of compressed sampling can lead to significant re- 
ductions in the computation of non-linear 2-D spectra. 
In such cases only a small number of data points need 
to be determined while the rest is constructed via matrix 
reconstruction techniques. Our findings are more gen- 
eral however and can provide considerable computational 
savings in a wide variety of computational problems in 
physics where data are represented by low-rank matrices. 

Acknowledgements - We acknowledge discussions with 
T. Baumgratz, F. Caycedo-Soler, M. Cramer and S.F. 
Huelga. This work was supported by an Alexander 
von Humboldt Professorship, the EU Integrated Project 
QESSENCE, the BMBF Verbundprojekt QuOReP, the 
Ministerio de Ciencia e Innovation Project No. FIS2009- 
13483-C02-02 and the Fundacion Seneca Project No. 



5 



11920/PI/09-j. 



[1] E.J. Candes and M.B. Wakin, An introduction to com- 
pressive Sampling, IEEE Signal Processing Magazine 25 
2008, 21-30 

[2] E. Candes, J. Romberg and T. Tao, Robust uncertainty 
principles: Exact signal reconstruction from highly in- 
complete frequency information, IEEE Trans. Inform. 
Theory. 2006, 52, 489 - 509. 

[3] E. Candes and T. Tao, Near optimal signal recovery from 
random projections: Universal encoding strategies, IEEE 
Trans. Inform. Theory. 2006, 52, 5406 - 5425. 

[4] D. Donoho, Compressed Sensing, IEEE Trans. Inform. 
Theory. 2006, 52, 1289 - 1 306. 

[5] See http://dsp.rice.edu/cs for a selection of review arti- 
cles on compressive sensing. 

[6] A. Ishizaki and G.R. Fleming, Unified treatment of quan- 
tum coherent and incoherent hopping dynamics in elec- 
tronic energy transfer: Reduced hierarchy equation ap- 
proach, J. Chem. Phys. 2009, 130, 234111 (10 pages). 

[7] J. Prior, A.W. Chin, S.F. Huelga and M.B. Plenio, Ef- 
ficient simulation of strong system-environment interac- 
tions, Phys. Rev. Lett. 2010, 105, 050404 (4 pages). 

[8] A.W. Chin, A. Rivas, S.F. Huelga and M.B. Plenio, Ex- 
act mapping between system-reservoir quantum models 
and semi-infinite discrete chains using orthogonal poly- 
nomials, J. Math. Phys. 2010, 51, 092109 (20 pages). 

[9] Ch. Kreisbeck and T. Kramer, Long-Lived Electronic 
Coherence in Dissipative Exciton-Dynamics of Light- 
Harvesting Complexes, E-print \arXiv: 1203. lj85\ 2Q12 (7 
pages). 

[10] X. Andrade, J. N. Sanders, A. Aspuru-Guzik, Applica- 
tion of compressed sensing to the simulation of atomic 
systems, E-print arXiv:1205.6485 2012 (7 pages). 

[11] J. N. Sanders, S. Mostame, S. K. Saikin, X. Andrade, J. 
R. Widom, A. H. Marcus, A. Aspuru-Guzik, Compressed 
sensing for multidimensional electronic spectroscopy ex- 
periments, E-print \arXiv:1207. 376b] 2012 (5 pages). 

[12] E.J. Candes and Y. Plan, Matrix completion with noise, 
Proc. IEEE 2010, 98, 925 - 936. 

[13] E. Candes and B. Recht, Exact Matrix Completion via 
Convex Optimization, Found. Comp. Math. 2009, 9, 717 
- 772. 

[14] D. Gross, Recovering Low- Rank Matrices From Few Co- 
efficients In Any Basis, IEEE Trans. Inf. Theo. 2011, 57, 
1548 - 1566. 

[15] J.-F. Cai, E.J. Candes and Z. Shen, A singular value 
thresholding algorithm for matrix completion, SIAM J. 



Opt. 2010, 20, 1956 - 1982. 

[16] M. Cramer, M.B. Plenio, S.T. Flammia, R. Somma, D. 
Gross, S.D. Bartlett, O. Landon-Cardinal, D. Poulin, 
Y.-K. Liu. Efficient quantum state tomography, Nature 
Comm. 2010, 1, 150. 

[17] S. Ma, D. Goldfarb, and L. Chen, Fixed point and Breg- 
man iterative methods for matrix rank minimization, J. 
Math. Prog. A and B 2011, 128, 321-353. 

[18] S. Becker, J. Bobin and E.J. Candes, NESTA: A Fast 
and Accurate First-Order Methods for Sparse Recovery, 
SIAM J. Imag. Set. 2011, 4, 1-39. 

[19] S. Boyd and L. Vandenberghe, Convex Optimization, 
Cambridge University Press 2004. 

[20] S. Mukamel, Principles of nonlinear optical spectroscopy, 
Oxford University Press 1995. 

[21] M. Cho, Two dimensional optical spectroscopy, CRC 
press, Taylor and Francis group 2009. 

[22] G. Engel, T. Calhoun, E. Read, T. Ahn, T. Mancal, Y. 
Cheng, R. Blankenship and G.R. Fleming, Evidence for 
wavelike energy transfer trough quantum coherence in 
photosynthetic systems, Nature 2007, 446, 782-786. 

[23] G. Panitchayangkoon, D. Hayes, K. Fransted, J. Caram, 
E. Harel, J. Wen, R. Blankenship, and G. Engel, Lon- 
lived quantum coherence in photosynthetic complexes at 
physiological temperature, Proceedings of the National 
Academy of Sciences 2010, 107, 12766-12770. 

[24] D. Hayes, G. Panitchayangkoon, K. A. Fransted, J. R. 
Caram, J. Wen, K. F. Freed, and G. S. Engel, Dynam- 
ics of electronic dephasing in the Fenna-Matthews-Olson 
complex, New J. Phys. 2010, 12, 065042 (12 pages). 

[25] F. Caycedo-Soler, A.W. Chin, J. Almeida, S.F. Huelga, 
and M.B. Plenio, The nature of the low energy band 
of the Fenna-Matthews-Olson complex: Vibronic signa- 
tures, J. Chem. Phys. 2012, 136, 155102 (13 pages). 

[26] B. Hein, C. Kreisbeck, T. Kramer, M. Rodriguez, Mod- 
elling of Oscillations in Two-Dimensional Echo-Spectra 
of the Fenna-Matthews-Olson Complex, New Journal of 
Physics 2012, 14, 023018 (20 pages). 

[27] F. Caruso, S. K. Saikin, E. Solano, S. F. Huelga, A. 
Aspuru-Guzik and M. B. Plenio, Probing biological light- 
harvesting phenomena by optical cavities, Phys. Rev. B, 
2012, 85, 125424 (10 pages). 

[28] D. E. Tronrud, J. Wen, L. Gay and R.E. Blankenship, 
The structural basis for the difference in absorbance spec- 
tra for the FMO antenna protein from various green sul- 
fur bacteria, Photosynth. Res. 2009, 100, 79-87, PDB 
ID:3EOJ 

[29] This time scale can be improved considerably by mak- 
ing use of the sparsity of the matrices and algorithmic 
improvements [TTl IT8] . 



